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Abstract — Traditional boundary integral methods suffer from 
the singularity of Green's kernels. The paper develops, for a 
model problem of 2D scattering as an illustrative example, 
singularity-free boundary difference equations. Instead of con- 
verting Maxwell's system into an integral boundary form first 
and discretizing second, here the differential equations are first 
discretized on a regular grid and then converted to boundary 
difference equations. The procedure involves nonsingular Green's 
functions on a lattice rather than their singular continuous 
counterparts. Numerical examples demonstrate the effectiveness, 
accuracy and convergence of the method. It can be generalized 
to 3D problems and to other classes of linear problems, including 
acoustics and elasticity. 

Index Terms — Scattering, diffraction, difference equations, 
boundary difference equations, boundary integral equations, 
boundary element methods, flexible local approximation, Green 
functions, discrete transforms. 

I. Introduction 

Boundary equation methods have a long history, with prac- 
tical applications dating back to the 1960s. An interesting 
historical account given by Cheng & Cheng [1] includes the 
work on wave scattering and radiation in 1962-1967 by Fried- 
man & Shaw, Chen & Schweikert, Banaugh & Goldsmith, 
Mitzner, and others [2]-[8]. In eletromagnetics, boundary 
integral techniques became very popular due primarily to 
Harrington's work published in 1967-68 [9], [10] (see also 
[HH15]). 

In traditional boundary integral methods, linear boundary 
value problems of field analysis are transformed into integral 
equations with respect to equivalent sources residing on the 
boundaries. In the simplest example of capacitance calculation 
[9], [10], the distributed charge density on conducting plates 
becomes the principal unknown variable. By equating the 
Coulomb potential of that charge to the given potential of 
the conductors, one obtains an integral equation. It can then 
be discretized using variational techniques (moment methods), 
collocation and Galerkin methods being particular cases of 
those. 

As all numerical methods, boundary integral techniques 
do carry some trade-offs. Their key advantage is the lower 
dimensionality of the problem: 3D analysis is reduced to 
equivalent sources on 2D boundaries and 2D analysis - to 
ID contours. Another advantage is a natural treatment of 
unbounded problems (e.g. wave scattering and radiation), with- 
out the artificial domain truncation unavoidable in differential 
methods such as finite difference (FD) schemes and the Finite 
Element Method (FEM). 



Integral equation methods have, in general, two major disad- 
vantages. First, the matrices of the discrete system are almost 
always full. This is due to the fact that a source at any point 
on the boundary contributes to the field at all other points. In 
contrast, FD and FE matrices are sparse, with very efficient 
system solvers available (iterative: multilevel methods, incom- 
plete factorization and other effective preconditioners; direct: 
minimum degree, nested dissection and others; see e.g. [16] 
and references there). Cases where Green's functions decay 
rapidly in space, giving rise to quasi-sparse integral equations, 
are exceptional (e.g. periodic structures in the electromagnetic 
band gap regime [17]). 

Another disadvantage is that the integral kernels in field 
analysis are singular. At the surface points, the kernel sin- 
gularity can usually be handled analytically, and the fields 
remain bounded as long as the surfaces are smooth. However, 
for points in the vicinity of the surface, the evaluation of the 
integral is problematic, as analytical expressions are usually 
unavailable and numerical quadratures require extreme care. 
The same is true for two adjacent surfaces with a narrow gap 
in between. 

Significant progress in Fast Multipole Methods (FMM) [18], 
[19], [20], [21], [22] has helped to alleviate the first disadvan- 
tage of boundary methods. FMM accelerates the computation 
of fields due to distributed sources - or equivalently, matrix- 
vector multiplications for the dense system matrices. 

The second disadvantage is more difficult to overcome. 
Singular kernels are inherent in boundary integral methods 
because the fields of point sources are unbounded. However, 
a drastic change in the computational procedure leads to a 
singularity-free method; this is accomplished by reversing the 
sequence of stages in the boundary techniques. The standard 
sequence is 

Differential formulation =>• Boundary 
integral formulation =>• Discretization 

The alternative sequence is 

Differential formulation =^> 
Discretization =>• Boundary difference 
problem 

Discretization of the differential problem is performed on 
a regular grid and yields an FD scheme. This scheme is 
converted - as explained in the remainder of the paper - to a 
boundary problem that involves discrete fundamental solutions 



(Green's functions) on the grid. Discrete Green's functions, 
unlike their continuous counterparts, are always nonsingular. 

This general idea is not new. In fact, there are two related 
but independently developed methodologies for boundary dif- 
ference equations. The first one, put forward and thoroughly 
studied by Ryaben'kii, Reznik, Tsynkov and others [34], [35], 
[36], [37], is known as the method of difference potentials and 
can be viewed as a discrete analog of the Calderon projection 
operators in functional analysis [34]. 

The second methodology, called boundary algebraic equa- 
tions by Martinsson & Rodin [23], is at least 50 years old 
(Saltzer [25]) and is a discrete analog of first- or second-order 
Fredholm boundary integral equations for potential problems 
[23]. 

In comparison with [23], the method of this paper has 
several novel features. First, the paper deals - to my knowl- 
edge, for the first time - with boundary difference equations 
for electromagnetic wave scattering. In [23], a simple model 
problem is considered: the Laplace equation (e.g. electrostatics 
or heat transfer) in a homogeneous domain with Dirichlet 
boundary conditions; the focus of [23] is on the mathematical 
analysis of the respective boundary difference operators, their 
spectral properties and the appropriate iterative solvers. 

One key distinction between the methodology of Ryaben'kii 
[34] and this paper's is in the choice of the main unknown: 
the boundary field / potential (Ryaben'kii) vs discrete sources 
on the boundary (the present paper). The treatment via sources 
parallels that of the continuous boundary integral method [10], 
[11], [12], [14], [15] and should therefore be intuitive to 
applied scientists and practitioners. Further analysis and com- 
parison of these methodologies will be presented elsewhere. 

An additional novelty of this paper is the use of high- 
order Flexible Local Approximation MEthods (FLAME, see 
Section IV- A) in the context of boundary difference equations. 
Also, this is the first application of FLAME to a 2D boundary 
of a generic shape; this is done by approximating this boundary 
locally by its osculating circle at any given point. 

II. Boundary Difference Equations for a Model 
Problem 

A. Formulation and Setup 

To fix ideas and explore the potential of the proposed 
approach, let us consider the classical 2D case of electro- 
magnetic wave scattering as a model problem. It should be 
emphasized from the outset that the method has a much 
broader range of applicability; possible generalizations are 
discussed in Section V-C. 

Consider a plane wave impinging from the air on a dielectric 
cylinder (Fig. 1) with a given dielectric permittivity e cy \. The 
cross-section of the cylinder could be arbitrary, but for the 
sake of simplicity we shall assume that its surface is smooth 
(no edges or corners). 

For definiteness, let us focus on the E-mode (TM- or s- 
mode) governed by the familiar equation for the electric field 
E with a single z-component: 

V 2 E(x,y) + k 2 (x,y)E(x,y) = 0, k 2 = uj 2 jie (1) 



E 




Fig. 1. Setup of the scattering problem for the E-mode. 

where the standard notation for the angular frequency u, the 
magnetic permeability /i, the dielectric permittivity e and the 
wavenumber k is used (k is equal to k cy \ inside the scatterer 
and to k out outside). Equation (1) should be supplemented by 
the standard radiation boundary conditions for the scattered 
field E s — E — E- inc at infinity. The incident field is a plane 
wave 

£inc = £ exp(-jk out • r), r=(x,y) (2) 

where the ex.p(+jujt) convention for complex phasors is 
implied. 

In a departure from the boundary integral methodology, we 
now proceed, prior to formulating a problem on the boundary 
of the scatterer, to FD discretization. To this end, let us 
introduce an infinite lattice with a grid size h, for simplicity 
the same in both x and y directions. Although infinite lattices 
are not a very common computational tool, they were already 
featured prominently in Martinsson & Rodin's work [23], [39] 
as well as in the much earlier report by Saltzer [25]. The actual 
computation, clearly, never involves an infinite amount of data 
on the lattice; in fact, the unknowns are ultimately confined 
only to the boundaries. 

As an auxiliary device, we need to consider the wave 
equation (1) in the homogeneous space with a constant generic 
parameter k. Various FD discretizations of this equation are 
available; see e.g. Harari's review [24] for further information 
and references. Here we settle for the simplest five-point 
scheme 

C(h,k)E = E(m x — l,m y )+E(m x + l,m y )—4:E(rn x ,rn y ) 

+E(m x ,m y -l) + E(m x ,m y + l) + k 2 h 2 E(m x ,m y ) = 

(3) 

where E(m Xl m y ) is the field value at a grid point charac- 
terized by an integer double index m = (m x ,m y ) G Z 2 . 
As reflected in the notation, the coefficients of the difference 
operator C depend on the mesh size and on the wavenumber; 
this may not be explicitly indicated if there is no possibility 
of confusion. 

Associated with C is its Green's function g(m XJ m y ) defined 
as the solution of 

Cg = S (4) 
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Fig. 2. Discrete boundary 7 with 196 nodes. Squares: 7i n ; circles: 7out- 

with the boundary condition 

g(m x ,m y ;k) — » G(m x h,m y h;k) as (ra^, m :y ) — >> 00 (5) 

Here J is the discrete delta-function (equal to one at the origin 

(2) 

and zero elsewhere) and G(r; fc) = ffi ; (/cr) is the continuous 

(2) 

Green function, Hq ' being the Hankel function. 

Without getting into the mathematical theory of lattice 
Green functions (see [39], [20] and Section III), let us note 
some features critical for our purposes: 

• In contrast with its continuous counterpart G, the discrete 
Green function (4) is bounded everywhere, including the 
origin. 

• The discrete Green function differs significantly from the 
continuous one only within a spatial window of several 
grid layers around the origin. Therefore only a relatively 
small amount of information needs to be stored - namely, 
the values of the Green function within that window. This 
data can be precomputed for any given value of kh and 
for each linear medium in a given problem. 

The discrete boundary of the scatterer can be defined in a 
natural way follwoing Ryaben'kii [34]. Each grid node m with 
discrete coordinates (m x ,m y ) has four immediate neighbors 
from the respective five-point stencil in the difference scheme 
(3). If node m lies inside the scatterer but some of his 
neighbors are outside, this grid node will be said to belong 
to the discrete inner boundary 7i n . Likewise, if the central 
node m of the stencil lies outside the scatterer but at least one 
of its neighbors is inside, this node is said to belong to the 
discrete outer boundary 7 out . The complete discrete boundary 
consists of two layers: 7 = 7 in U 7 ou t> Fig- 2. (For larger 
multipoint stencils, the discrete boundary can be composed 
of several layers.) The number of nodes on the boundary is 
n 7 = n 7) i n + n 7)OU t- These nodes can be referred to by pairs 
of indexes (m x , m y ) or, alternatively, by some global numbers 
from 1 to n 7 . The order of this numbering makes no principal 
difference but may slightly affect the practical implementation 
of the method. 

B. Boundary Sources 

The critical step is to express the lattice-based field in terms 
of fictitious discrete sources / that are nonzero only on the 



discrete boundary 7. For the inner boundary, 

E(m) = [/*0(v;fccyi)](m), m = (m x ,m y ) e 7in (6) 
For the outer boundary, 

E(m) = E inc (m) + [/ * #(•, fc out )](m), m £ 7 out (7) 

The discrete convolution in the equations above is defined in 
the usual way as 

[f*g](n) = ^/(m) 5 (n-m) (8) 

mG7 

Note that for the nodes on each side of the boundary the field 
is described via the respective discrete Green function. 

The auxiliary sources / need not have a direct physical 
interpretation, although ultimately they are indirectly related 
to the equivalent electric and magnetic surface currents of 
traditional boundary integral methods [14], [13]. However, the 
fields derived from these sources are physical. The convo- 
lutions in equations (6), (7) can be interpreted as (discrete) 
scattered fields. 

It can be shown that any discrete field satisfying the FD 
wave equation on the lattice can indeed be expressed via 
convolution with some fictitious boundary sources as stipu- 
lated above, except possibly for special resonance cases (see 
Appendix). 

C. Boundary Difference Equations 

By construction, the electric fields defined by (6), (7) satisfy 
the respective wave equation on each side of the boundary. 
What remains to be done then is to impose the boundary 
conditions; this will lead to a system of equations from which 
the sources can be found. 

To this end, one may use another difference scheme, <S, 
that approximates the boundary conditions; we shall call it a 
"boundary test scheme." The simplest example is the five-point 
scheme 

S 5 E E(m x — l,m y ) + E(m x + l,m y ) —4E(m x ,m y ) 

+E(m Xj m y -l)+E(m x ,m y +l)+k 2 (m)h 2 E(m) = (9) 

In this second-order scheme, the value of k is taken at the 
midpoint of the stencil. A more accurate alternative is the 
nine-point FLAME (Flexible Local Approximation MEthod) 
proposed in [26], [27], [16]). Both types of schemes are 
used in the numerical examples of Section IV. The FLAME 
coefficients <Sg are computed as the nullspace of a matrix 
comprising the nodal values of a set of basis functions on 
the stencil (Section IV-A and [26], [27], [16]). 

Applying a given boundary test scheme S on 7 to fields (6), 
(7), one obtains a system of boundary-difference equations of 
the form 

5 ( ™ ) [/* 5 (-,-;fc(n)](n) = -S^E inC:h ; n=(n x ,n y )ej 

(10) 

where E- inCt h is the discrete version of the incident field (i.e. 
its values on the lattice). The superscript (n) indicates that 
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different schemes with different coefficients can in principle 
be used over different stencils. 

More explicitly, denoting the coefficients of the boundary 
test scheme with (where index a runs over nodes 
a over the grid stencil centered at node n), one can write the 
boundary equation (10) as 

Af = q (11) 

where A is an n 7 x n 7 matrix with the entries 

a 

and 

a - -V^E M 

( in — / j ol ^inc,^ 

OL 

The meaning of the terms above is as follows: 

• n, m are the global numbers 1 (1 < n,m < n 7 ) of 
nodes n = (n x ,n y ) and m = (m x ,m y ) on the discrete 
boundary 7. 

• are the coefficients of the boundary test scheme 
corresponding to node n. (In principle, different schemes 
could be used at different nodes. One may even envision 
an adaptive procedure where the order of the scheme will 
vary in accordance with local accuracy estimates.) 

• k = k cy \ if node n is on the inner boundary 7i n and 
k = k out if it is on the outer boundary 7 out . 

• ^inch * s me vame °f the incident field at node a of 
stencil n. 

We shall call the numerical procedure leading to (11) the 
boundary difference method (BDM). 

III. The Lattice Green Function 

As evident from the description of the BDM, lattice Green's 
functions play a central role in it and must be computed ac- 
curately. There are at least two general ways to do so: Fourier 
analysis and finite difference solutions. A detailed exposition 
for the Laplace equation has been given by Martinsson & 
Rodin [20], [39], [23]. Similar ideas can be immediately ap- 
plied to the wave equation as well, although a more elaborate 
analysis would be desirable in the future. 

Applying Fourier transform T (discrete physical space — )> 
continuous reciprocal space) to the difference equation (4) with 
the five-point operator £ (3), one obtains 

T{Cg) = (exp(jft x ) + exp(-j^) + exp(j^) + exp(-j^) 

-4 + k 2 h 2 )T{g} = 1 

where n x , n y are the Fourier parameters in the square [— 7r, 7r] 2 . 

The inverse Fourier transform may then serve as a staring 
point for an asymptotic analysis similar to Martinsson's [39], 
[20] and for practical computation of Green's function g. 
However, this Fourier analysis is quite involved and must be 
performed with great care, especially in 2D where Green's 
functions decay slowly and regularization of Fourier integrals 
is necessary [39], [20]. For the purposes of this paper, a 

! Not to be confused with the Euclidean lengths of n = (n x ,n y ) and 
m = (m x ,m y ); these lengths are irrelevant and never appear in our analysis. 
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Fig. 3. Re(#) for k = 1, h = 1/7, M = 50. Note that the discrete green 
function is nonsingular everywhere; in fact, its magnitude in this example is 
quite moderate. 
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Fig. 4. Re(#) for k = 2, h = 1/7, M = 50. 

more straightforward route is sufficient. The finite difference 
problem (4) for the Green function can be solved directly, with 
the boundary condition (5) imposed on the boundary of a large 
enough square [— M, M] 2 . This can be done efficiently with 
fast Fourier transforms over the square, but the computational 
cost in 2D is so moderate that any other reasonable solver 
can be applied. Obviously, one can also take advantage of the 
symmetries to reduce the size of the computational problem. 

The following plots illustrate the behavior of the lattice 
Green function and its computation. All of the plots were 
generated for the grid size h = 1/7 as an example. Surface 
plots of the real part of g for wavenumbers k = 1 and k = 2 
are shown in Figs. 3 and 4, respectively; Green's function was 
computed in the spatial window [— M, M] 2 with M = 50. 

Fig. 5 demonstrates that the size M of the window need not 
be too large. Indeed, lattice Green's functions for M = 50 and 
M = 100 are quite close. The numerical experiments reported 
in the following section were performed with M = 50. 
Even assuming an overkill value M = 100 and 10 different 
materials in a given practical problem, one ends up with 
less than 1 MB of data to be stored. In 3D, if one takes 
advantage of the symmetries of g, the memory requirements 
are still reasonable, even for vector fields and dyadic Green's 
functions, except for problems where the number of different 
materials is unusually large. 

IV. Numerical Simulation 

A. FLAME 

The theory, implementation and various applications of 
FLAME have been discussed in a number of previous pub- 
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Fig. 5. Lattice Green's function g(m x ,m y ) vs m x for kh = 1/7, m y = 0. 
The results for two different values of M, M = 50 and M = 100, are close. 



lications [26], [27], [16], [32], [33], [29], [31], and therefore 
only a brief summary is given here. 

FLAME replaces the usual Taylor expansions, the key tool 
of standard finite-difference analysis, with much more accurate 
local (quasi-)analytical approximations of the solution. Such 
approximations can be obtained, for example, via cylindrical or 
spherical harmonics, plane waves, etc. Since the local behavior 
of the field is "built into" the difference scheme, the accuracy 
often improves dramatically. FLAME has already been applied 
to the simulation of colloidal and plasmonic particles [26], 
[27], [16], negative-index materials [30], the computation of 
Bloch bands of photonic crystals [28], including complex 
bands for plasmonic systems and other dispersive media [28], 
[29]. 

For the model problems in this paper, local analytical bases 
for FLAME are available via Bessel / Hankel functions. More 
specifically, in the vicinity of a dielectric cylinder with a cir- 
cular cross-section centered for convenience at the origin of a 
polar coordinate system (r, </>), these approximating functions 
- the FLAME basis ^ - are [26], [27], [16] 



a^Ji(k cy \r) exp(j7(/>), r < r cy i 



(2) 

r ' Jl(k^ T r) + H\ H/CoutO] ex P(j#)> r > r cyl 

(2) 



where J\ is the Bessel function of order H\ ' is the Hankel 

(i) (i) 

function of the second kind, and a\ , c\ ' are coefficients to 
be determined. These coefficients are found via the standard 
conditions on the boundary of the cylinder [26], [16], [28]. 
Index i runs over all grid stencils where the FLAME scheme 
is generated, while index a runs over all basis functions in a 
given stencil i. 

In this paper, the 9-point (3x3) stencil with a grid size h is 
used and 1 < a < 8. The eight basis functions ip are obtained 
by retaining the monopole harmonic (/ = 0), two harmonics 
of orders \l\ = 1,2,3 (i.e. dipole, quadrupole and octupole), 
and one of the harmonics of order |2| =4. This set of basis 
functions produces a nine-point scheme as the null vector of 
the respective matrix of nodal values [26], [27]. 

For the test problem with an elliptical cylinder (Sec- 
tion IV-C), it is still possible to use the same Bessel-Hankel 




Fig. 6. Color plot of the E field for a circular cylinder with e cy \ = 4. BDM 
with n 7 = 460. 



basis functions in FLAME. Toward this end, a piece of the 
ellipse straddled by a given grid stencil is approximated by 
its osculating circle (with the radius equal to the radius of 
the curvature of the ellipse at a given point on its boundary). 
While this approach for constructing FLAME bases is fairly 
straightforward, it has never been used previously. (In the past, 
the primary motivation was to apply FLAME on very coarse 
grids that carry almost no information about the shape of the 
boundary [26], [27], [16].) 

For the ellipse, the osculating circle can easily be found 
analytically; for more complicated boundaries, the curvature 
could approximately be evaluated numerically - for example, 
as the best local fit to a piece of the discrete boundary 7. 
In yet more complex cases - especially in 3D where there 
are two radii of curvature - one could use piecewise-planar 
approximations and Fresnel-formula FLAME bases [30]. 

B. Circular cylinder 

For verification, let us first consider a circular scattering 
cylinder, as in this case a well-known analytical solution 
via cylindrical harmonics exists. In all numerical experiments 
below, the E mode was considered. The wavenumber for the 
incident wave was normalized at k out = 1; the wavenumber 
for the scatterer was taken as k cy \ = 2 (i.e. e cy i = 4). The 
incident plane wave propagates in the positive x direction. 
The color plot of the electric field in BDM with n 7 = 460 is 
shown in Fig. 6. 

The numerical error as a function of the BDM grid size 
h is plotted in Fig. 7 for two boundary test schemes S\ the 
standard five-point scheme and the nine-point FLAME (these 
schemes were briefly described in the previous sections). The 
dashed line in the figure serves only as a visual aid indicating 
the second order convergence of the method for both schemes. 
Not surprisingly, the numerical error for FLAME is about an 
order of magnitude lower than that of the five-point scheme. 
However, the order of convergence is still limited by the 
second-order five-point difference scheme used to compute 
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Fig. 7. The relative error in BDM as a function of grid size. Discrete boundary 
7 with 196 nodes. Quadratic convergence, commensurate with the order of 
the scheme for the lattice Green function, is observed. The FLAME results 
are about an order of magnitude more accurate than for the standard five-point 
scheme, even around h ~ 0.05 where the FLAME data points exhibit some 
scatter. 

the discrete Green function (Section III). The relative error 
was calculated as ||£bdm - inexact ||/ Inexact ||, where E B dm 
and inexact are the numerical and the quasi-exact solutions on 
the grid, respectively; the norms are Euclidean. The quasi- 
exact solution was computed via the standard expansion into 
cylindrical harmonics up to order 50. 

C. Elliptical Cylinder 

The simulations have been repeated for an elliptical cylin- 
der, with the same physical parameters as above, and with the 
ratio of the axes 1.5 : 1. Fig. 8 is a color plot of the real 
part of the electric field obtained with BDM, 9-point FLAME 
scheme, discrete boundary 7 with 196 nodes. FLAME was 
generated as described at the end of Section IV-A: by locally 
approximating a piece of the elliptic boundary with a circle 
and using the respective Bessel / Hankel bases. 

Fig. 9 demonstrates that the field distributions obtained with 
different methods are in a very good agreement. Plotted in the 
figure is the real part of the electric field vs. x (at y = 0) 
and vs y (at x = 0). The imaginary parts are not plotted, but 
agree with the theory equally well. Nine-point FLAME and 
standard five-point schemes were applied on several grids with 
sizes h — r cy i/(n + \)\ results for n = 10 and n = 20 are 
shown. 

V. Discussion and Conclusion 

A. Summary 

The boundary difference method described and imple- 
mented in this paper for wave scattering avoids the sin- 
gularities inherent in traditional boundary integral methods. 
This is accomplished by reversing the sequence of stages 
in the procedure. Traditionally, the differential equations are 
first reduced to boundary integrals with respect to equivalent 




x 



Fig. 8. Rq(E) obtained with BDM, 9-point FLAME scheme. Discrete 
boundary 7 with 196 nodes. 




n=20 

coordinate 

Fig. 9. The numerical results for different cases are seen to be in a very 
good agreement. The real part of the electric field is plotted; the agreement 
for the imaginary part is similar. Discrete boundary 7 with 196 nodes. 



sources on the boundary and then discretized; the kernels 
of the underlying integral equations are singular due to the 
infinite self-fields of concentrated sources. 

In BDM, the differential problem is first discretized on a 
regular grid to obtain a finite-difference approximation that is 
then reduced to a boundary difference equation with respect to 
auxiliary sources on the discrete boundary. The field of these 
sources can be expressed by convolution with the discrete 
Green function that, unlike its continuous counterpart, is finite 
at all points. Thus no singularities ever arise. 

Technically, the underlying grid is infinite. The computa- 
tional procedure, however, involves only the boundary nodes 
of the grid and a finite spatial window where the discrete Green 
function is precomputed, which can be done once and for all 
for a given set of parameters. 

The validity of BDM has been demonstrated using 2D 
scattering from dielectric cylinders as a model problem. Con- 
vergence of the method as a function of the grid size has 
been established and is commensurate with the order of finite- 
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difference schemes used. 

B. Trade-offs 

Since the proposed approach has common features with the 
traditional integral equation methods, some of the usual trade- 
offs between differential and integral techniques [38] apply. 
The differential methods lead to sparse matrices, whereas 
the boundary methods produce dense ones. This drawback 
can be partly alleviated via fast multipole acceleration [18], 
[20], [21], [23], [22]. Its use in conjunction with BDM is 
relatively straightforward. Indeed, FMM relies on a recursive 
splitting of the solution into near- and far-field components. 
The far field in BDM is essentially the same as in the 
continuous problem, by construction of the discrete Green 
function; see (5). It is only in the near field that discrete and 
continuous Green functions may differ significantly, but this 
makes little difference in FMM algorithms because the near- 
field contribution is computed directly. 

As already emphasized, BDM completely dispenses with 
singular integral kernels, an inherent drawback of integral 
methods. The price to pay for that is the need to precompute 
discrete Green functions. In practice, this price can be expected 
to be modest, because the number of different materials in 
any given problem is limited and the computation involves 
a relatively small number of grid layers around Green's 
point source. In any event, this computational overhead is 
independent of the size of the problem being solved. 

For unbounded problems, differential methods such as FEM 
and FD require artificial domain truncation with absorbing 
boundary conditions or matched layers. No such truncation is 
needed in boundary methods. At the same time, differential 
methods are generally better suited for nonlinear problems 
that call for volume discretization, in which case the boundary 
methods usually lose their effectiveness. 

The key source of numerical errors in traditional BEM is 
approximation of singular integrals (typically, by piecewise- 
polynomial functions of low order, including piecewise- 
constant approximations in the simplest case). In BDM, the 
error is due to the finite-difference approximation of the 
boundary conditions and of the lattice Green function. If 
the order of these approximations is increased, the overall 
numerical error of the method can be reduced accordingly. 

C. Generalizations and Future Directions 

Boundary difference schemes developed in this paper lend 
themselves to generalization in quite a natural way. Unlike 
traditional boundary integral methods, BDM is automatic, 
in the sense that it does not require the suitable sets of 
equivalent boundary sources (electric or magnetic surface 
currents, surface charges, etc.) and the respective equations 
to be worked out in advance. Instead, one introduces discrete 
boundary sources that need not even have a specific physical 
meaning; but once computed, they can be used to find physical 
fields by convolution with Green's functions on the lattice. In 
particular, the i^-mode (TE- or p-mode) of electromagnetic 
wave scattering is treated in BDM exactly the same way as 
the £?-mode. (As a side note, for the i^-mode the classical 



five-point control volume scheme would only be of order one 
at the boundary, but that is a feature of that scheme, not of 
BDM as a whole.) 

Further, extension to 3D vector problems is also conceptu- 
ally straightforward, although clearly algorithmic challenges 
do arise. This line of research is currently being pursued. 

The boundary difference method does in general require a 
spatially uniform grid. Although this grid is "virtual," in the 
sense that the actual computation involves only the nodes on 
the discrete boundary and not the volume nodes, the uniformity 
of the grid may still be a limiting factor in some problems. 
However, if several scatterers are present and well- separated 
(in practice, by at least a few grid layers), then each of them 
may be meshed separately. Indeed, in that case the interactions 
between different scatterers are numerically in the far field, 
where the continuous Green function can be used as a good 
proxy for the discrete one. 

Finally, the method is not limited to electromagnetics and 
can be extended to other classes of linear problems, including 
acoustics and elasticity. It may even be applied to micro-, 
nano- and molecular- scale models on a discrete lattice (e.g. 
Haq et al. [40]), when continuous equations may not even be 
available. 

Appendix: Representation of the Field via 
Discrete Sources 

Let us show that any discrete field on the boundary 7 can be 
represented via convolution of the Green functions with some 
auxiliary sources on the same boundary, except possibly for 
some special cases of interior resonance. More precisely, let 
E s h = Eh—Ei nc ^ be the scattered component of a lattice field 
Eh that satisfies the discretized wave equation both inside and 
outside the scatterer. Further, let E s h n represent the values of 
E s h on the discrete boundary 7. We intend to show that 

E shn (m) = [/*^(.,.;fc(m))](m) (12) 

for some source / on 7. 

The discrete convolution in (12) can be viewed as a linear 
operator that maps functions / in R n ^ to fields E s h n , also in 
R n ~i . It is then sufficient to demonstrate that this operator is 
nonsingular or, equivalently, that an identically zero field on 
the discrete boundary can be produced only by zero sources 
on that boundary. 

Let us thus assume that E s h n is identically zero. Then E s h 
must be zero, too, everywhere in the outside region. This is true 
because, by its construction as convolution with sources only 
on 7, this field satisfies the homogeneous difference equation 
in the outside region and also, by assumption, the Dirichlet 
conditions for it on 7 out are zero. Similar considerations hold 
for E s h in the inside region away from the interior resonance, 
as long as k cy \ is not an eigenvalue of the wave problem 
inside the scatterer, with zero Dirichlet conditions. Thus the 
convolution in (12) must be identically zero on the whole 
lattice, from which it immediately follows (e.g. via Fourier 
transforms) that / = 0. 
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